! Copyright (c) [2019] Bangjie Deng
! FChaseSolver is licensed under the Mulan PSL v1.
! You can use this software according to the terms and conditions of the Mulan PSL v1.
! You may obtain a copy of Mulan PSL v1 at:
! http://license.coscl.org.cn/MulanPSL
! THIS SOFTWARE IS PROVIDED ON AN "AS IS" BASIS, WITHOUT WARRANTIES OF ANY KIND, EITHER EXPRESS OR
! IMPLIED, INCLUDING BUT NOT LIMITED TO NON-INFRINGEMENT, MERCHANTABILITY OR FIT FOR A PARTICULAR
! PURPOSE.
! See the Mulan PSL v1 for more details.


submodule (FChaseSolver) FChaseSolverImpl
    implicit none;
contains
    module  subroutine fcs_idep(x, n, a, b, c, d)
        integer, parameter :: k = 4;
        integer, intent(in) :: n;
        real(kind = k), dimension(:), intent(in) :: a(:);
        real(kind = k), dimension(:), intent(in) :: b(:);
        real(kind = k), dimension(:), intent(in) :: c(:);
        real(kind = k), dimension(:), intent(in) :: d(:);
        real(kind = k), dimension(:), intent(inout) :: x;


        real, allocatable :: l(:);
        real, allocatable :: f(:);

        integer :: i;

        allocate(l(n));
        allocate(f(n));

        l(1) = b(1);
        do i = 2, n
            l(i) = b(i) - a(i - 1) * c(i - 1) / l (i - 1)
        end do

        f(1) = d(1) / l(1)
        do i = 2, n
            f(i) = (d(i) - a(i - 1) * f(i - 1)) / l(i)
        end do

        x(n) = f(n)
        do i = n - 1, 1, -1
            x(i) = f(i) - x(i + 1) * c(i) / l(i)
        end do

        deallocate(l);
        deallocate(f)

    end subroutine fcs_idep;

    module  subroutine fcd_idep(x, n, a, b, c, d)
        integer, parameter :: k = 8;
        integer, intent(in) :: n;
        real(kind = k), dimension(:), intent(in) :: a(:);
        real(kind = k), dimension(:), intent(in) :: b(:);
        real(kind = k), dimension(:), intent(in) :: c(:);
        real(kind = k), dimension(:), intent(in) :: d(:);
        real(kind = k), dimension(:), intent(inout) :: x;


        real, allocatable :: l(:);
        real, allocatable :: f(:);

        integer :: i;

        allocate(l(n));
        allocate(f(n));

        l(1) = b(1);
        do i = 2, n
            l(i) = b(i) - a(i - 1) * c(i - 1) / l (i - 1)
        end do

        f(1) = d(1) / l(1)
        do i = 2, n
            f(i) = (d(i) - a(i - 1) * f(i - 1)) / l(i)
        end do

        x(n) = f(n)
        do i = n - 1, 1, -1
            x(i) = f(i) - x(i + 1) * c(i) / l(i)
        end do

        deallocate(l);
        deallocate(f)

    end subroutine fcd_idep;

    module  subroutine fcc_idep(x, n, a, b, c, d)
        integer, parameter :: k = 4;
        integer, intent(in) :: n;
        complex(kind = k), dimension(:), intent(in) :: a(:);
        complex(kind = k), dimension(:), intent(in) :: b(:);
        complex(kind = k), dimension(:), intent(in) :: c(:);
        complex(kind = k), dimension(:), intent(in) :: d(:);
        complex(kind = k), dimension(:), intent(inout) :: x;



        real, allocatable :: l(:);
        real, allocatable :: f(:);

        integer :: i;

        allocate(l(n));
        allocate(f(n));

        l(1) = b(1);
        do i = 2, n
            l(i) = b(i) - a(i - 1) * c(i - 1) / l (i - 1)
        end do

        f(1) = d(1) / l(1)
        do i = 2, n
            f(i) = (d(i) - a(i - 1) * f(i - 1)) / l(i)
        end do

        x(n) = f(n)
        do i = n - 1, 1, -1
            x(i) = f(i) - x(i + 1) * c(i) / l(i)
        end do

        deallocate(l);
        deallocate(f)

    end subroutine fcc_idep;

    module  subroutine fcz_idep(x, n, a, b, c, d)
        integer, parameter :: k = 8;
        integer, intent(in) :: n;
        complex(kind = k), dimension(:), intent(in) :: a(:);
        complex(kind = k), dimension(:), intent(in) :: b(:);
        complex(kind = k), dimension(:), intent(in) :: c(:);
        complex(kind = k), dimension(:), intent(in) :: d(:);
        complex(kind = k), dimension(:), intent(inout) :: x;


        real, allocatable :: l(:);
        real, allocatable :: f(:);

        integer :: i;

        allocate(l(n));
        allocate(f(n));

        l(1) = b(1);
        do i = 2, n
            l(i) = b(i) - a(i - 1) * c(i - 1) / l (i - 1)
        end do

        f(1) = d(1) / l(1)
        do i = 2, n
            f(i) = (d(i) - a(i - 1) * f(i - 1)) / l(i)
        end do

        x(n) = f(n)
        do i = n - 1, 1, -1
            x(i) = f(i) - x(i + 1) * c(i) / l(i)
        end do

        deallocate(l);
        deallocate(f)

    end subroutine fcz_idep;

    subroutine ccs_idep(x, n, a, b, c, d) bind(c)
        use, intrinsic :: iso_c_binding;
        implicit none;
        integer(kind = C_INT), intent(in), value :: n;
        REAL (c_float), TARGET, INTENT(IN) :: a(n), c(n)
        REAL (c_float), TARGET, INTENT(IN) :: b(n), d(n)
        REAL (c_float), TARGET, INTENT(OUT) :: x(n)

        call fcs_idep(x, n, a, b, c, d);
    end subroutine

    subroutine ccd_idep(x, n, a, b, c, d) bind(c)
        use, intrinsic :: iso_c_binding;
        implicit none;
        integer(kind = C_INT), intent(in), value :: n;
        REAL (c_double), TARGET, INTENT(IN) :: a(n), c(n)
        REAL (c_double), TARGET, INTENT(IN) :: b(n), d(n)
        REAL (c_double), TARGET, INTENT(OUT) :: x(n)

        call fcd_idep(x, n, a, b, c, d);
    end subroutine

    subroutine ccc_idep(x, n, a, b, c, d) bind(c)
        use, intrinsic :: iso_c_binding;
        implicit none;
        integer(kind = C_INT), intent(in), value :: n;
        complex (c_float_complex), TARGET, INTENT(IN) :: a(n), c(n)
        complex (c_float_complex), TARGET, INTENT(IN) :: b(n), d(n)
        complex (c_float_complex), TARGET, INTENT(OUT) :: x(n)

        call fcc_idep(x, n, a, b, c, d);
    end subroutine

    subroutine ccz_idep(x, n, a, b, c, d) bind(c)
        use, intrinsic :: iso_c_binding;
        implicit none;
        integer(kind = C_INT), intent(in), value :: n;
        complex (c_double_complex), TARGET, INTENT(IN) :: a(n), c(n)
        complex (c_double_complex), TARGET, INTENT(IN) :: b(n), d(n)
        complex (c_double_complex), TARGET, INTENT(OUT) :: x(n)

        call fcz_idep(x, n, a, b, c, d);
    end subroutine
end submodule

